Subordinated diffusion and CTRW asymptotics 
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Anomalous transport is usually described either by models of continuous time random walks 
(CTRW) or, otherwise by fractional Fokker-Planck equations (FFPE). The asymptotic relation 
between properly scaled CTRW and fractional diffusion process has been worked out via various 
approaches widely discussed in literature. Here, we focus on a correspondence between CTRWs 
and time and space fractional diffusion equation stemming from two different methods aimed to 
accurately approximate anomalous diffusion processes. One of them is the Monte Carlo simula- 
tion of uncoupled CTRW with a Levy a-stable distribution of jumps in space and a one-parameter 
Mittag-Leffler distribution of waiting times. The other is based on a discretized form of a subor- 
dinated Langevin equation in which the physical time defined via the number of subsequent steps 
of motion is itself a random variable. Both approaches are tested for their numerical performance 
and verified with known analytical solutions for the Green function of a space-time fractional dif- 
fusion equation. The comparison demonstrates trade off between precision of constructed solutions 
and computational costs. The method based on the subordinated Langevin equation leads to a 
higher accuracy of results, while the CTRW framework with a Mittag-Leffler distribution of waiting 
times provides efficiently an approximate fundamental solution to the FFPE and converges to the 
probability density function of the subordinated process in a long-time limit. 

PACS numbers: 05.40.Fb, 05.10.Gg, 05.60.-k 02.50.-r, 02.50.Ey, 



In the ubiquity of complex systems ob- 
served in Nature, non-Gaussian fluctuations 
prevail and transport properties deviate from 
the standard theory of Brownian motion. 
Various facets of "anomalous diffusion" have 
been extensively studied over the past years 
by use of either continuous time random 
walks (CTRW) or fractional kinetic equa- 
tions. Recent development of simulation 
techniques based on Langevin equation in 
random subordinated time links uniquely 
stochastic trajectories of anomalous diffu- 
sion with probability density functions gov- 
erned by fractional Fokker-Planck equations 
and provides efficient tools to study CTRW 
and their asymptotics. Our work quantifies 
convergence of CTRW to anomalous diffu- 
sion process by using two schemes of model- 
ing: We compare simulations of the Montroll- 
Weiss-Scher model with a Mittag-Leffler dis- 
tribution of waiting times with subordinated 
Langevin dynamics. Efficiency of both meth- 
ods in reproducing analytical results is tested 
along with analysis of numerical accuracy and 
computational costs. 



I. INTRODUCTION 

The continuous time random walk (CTRW) concept, 
as introduced in pioneering works by Montroll-Weiss- 
Scher [l[ has been used over the decades for modeling dif- 
fusion processes on lattices, including anomalous trans- 
port. Commonly, the CTRW can be considered as a 
point process with reward [2|]. Within such a framework 
evolution of a random walker position is described by 
a sequence of independent, identically distributed ran- 
dom positive variables T„ which arc interpreted as wait- 
ing times between consecutive jumps of the walker. In 
the simplest uncoupled version of the CTRW scenario, 
the jumps of length x n and waiting times are indepen- 
dent from each other. Correspondingly, the position of a 
diffusive particle can be represented by a random sum of 
independent random variables 
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with N(t) staying for the counting random process which 
gives the (random) number of jumps up to a time t 
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A decoupled CTRW is Markovian only if the waiting 
time distribution ip(t) is exponential |3j, i.e. when the 
corresponding counting process N(t) is Poissonian. In 
this case, the CTRW process (QJ has time-homogeneous 
(stationary) increments, i.e. it is an infinite divisible 
compound Poisson process. So called dependent CTRW 



models (when the waiting times T n and the returns x 
are coupled) were first studied by Shlesinger et al. 
in order to place a physically realistic upper bound on 
particle velocities x n /T n . Furthermore, as discussed by 
Meerschaert et al. [5|, Q any coupled model at all for 
which in the long time scale (c — > oo) the convergence 
(c- 1 / a X[ ct ],c- 1 / , T[ ct ]) =*> (A(t),D(t)) is met will have 
one of two kinds of limits: Either the dependence dis- 
appears in the limit (because the waiting times and the 
jumps are asymptotically independent), or else the limit 
process is one of those derivable from the Shlesinger's 
CTRW model j4[. In such a case the counting process 
N(t), as directly related to the jump-times process T(n) 
by the relation 



N(t) ^ n <S> T(n) ^ t, 



(3) 



follows the inverse random time distribution [7H9J with 
the waiting time defined as the inter-jump time inter- 
val T(n) — T(n — 1). Put it differently, both processes 
T(n) and N(t) can be viewed as mutually inverse random 
functions leading to the equivalence 



Prob {T(n) < t} = Prob {N(t) ^ n} 



p n '(n'(t))dn'. 



(4) 



In the limit when n becomes a continuous parameter [7J- 
tSj, for which the variable T(n) is assumed to be dis- 
tributed according to a strictly asymmetric, one sided 
^-stable distribution L„i, the probability density (PDF) 
Pn(n(t)) can be obtained from the corresponding PDF of 
the random time T(n) 



Pn(n(t)) = -— [ l Vil (t';n)dt' 



(5) 



where lu,i(t;n) = dL v ^\(t\n)l dt. Accordingly, the long- 
time limit process X N u\ is then described by the proba- 
bility density [Ml, 



(6) 



p(x,t) = / p(x,T)p n (T,t)dT 
JO 

which possesses the scaling property 

p(x,t) =t- u / a p(xt- v/a ,l). 



(7) 



The time series Eq. (p} can be otherwise characterized 
by the probability of jumps x(t n +i) — x(t n ) and waiting 
times t n+ i — t n . In case that these both variables are 
statistically independent with the waiting time tp(t) and 
the jump length p(x) PDFs having tails (x — > oo, t — > oo) 
of the power-law type 



-( v +i) 



and 



ip(t) oc t 



p{x) oc |x|" (Q+1) , 



(8) 



(9) 



the Laplace-Fourier transform p(q, u) of the PDF p(x, t) 
takes the form H [11 



p{q,u) 



dsexp[-s(u v + \q\ a )]. (10) 



Here < v < 1 and < a < 2 are stability indices of the 
corresponding time and jump length distributions. 

The CTRW scheme as described by Eqs. © and © 
asymptotically leads to the situation when the evolu- 
tion of the probability density p(x, t) of finding a ran- 
dom walker at the position x after n steps performed 
up to time t can be described by a continuous integral 
p(x,t) = J Q dTp(x,T)p n (T,t) being the solution to the 
fractional Fokkcr-Planck equation 0, [lj, Ha | 
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with the initial condition p{x, 0) = 5(x). In the above 
equation oD t ~ v denotes the Riemann-Liouville fractional 
(time) derivative qD\~ v = ^D^ v defined by the rela- 
tion 
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and 



-^ stands for the Riesz-Weyl fractional (space) 

derivative with the Fourier transform T\ ai^if ] = 
— |fc| a J-"[/(a;)]. In the limit of a = 2 and v = 1, a Levy 
random walk ([]} is (asymptotically) equivalent to the 
standard (Markovian) Brownian motion, while for v = 1 
with a < 2 it corresponds to (Markovian) Levy flights. 

Both approaches are in use and have received atten- 
tion and application in a plethora of physical and bio- 
logical problems such as mass and charge transport in 
disordered systems, relaxation phenomena, front propa- 
gation in reaction-diffusion systems, transport in plasma, 
motion of organelles or epidemic spread [jj, [l7H26| . Suit- 
ability of the fractional dynamics by means of the cou- 
pled Langevin equations representing the subordination 
scheme has been also proved in analysis of correlation 
functions [27j |. 

Our paper is devoted to a detailed analysis of the 
CTRW diffusive limit obtained based on two different 
approaches. First relates to the role played by Mittag- 
Leffler functions in anomalous relaxation [lj, |28( . In par- 
ticular, we examine the uncoupled case of the CTRW sce- 
nario which includes implicitly the Mittag-Lcffler waiting 
time distribution function 
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1) 



being solution of the equation 

I>r*(T) - -*(r) 



(13) 



(14) 



for the so called survival function ^(r) = 1 — L ip(s)ds, 
^?(t) = E v {—t v ). The second choice is based on the 
CTRW modeling involving a subordination technique 

The paper is organized as follows: In Section |TT] 
some generic properties of both approaches are intro- 
duced. Section IIIII presents results of numerical analy- 
sis, in which asymptotic properties of the subordinated 
Langcvin equation are compared with the asymptotics of 
the CTRW scheme with ip(t) given by a one-parameter 
Mittag-Leffler distribution and p{x) represented by a 
symmetric Levy a-stable distribution. The supremacy of 
subordination technique is further demonstrated in Sec- 
tion IIVI which discusses numerical costs of computation 
and precision of both methods applied. The paper is 
closed with concluding remarks. 



II. RELATION BETWEEN CTRW, 

SUBORDINATION AND FRACTIONAL 

CALCULUS 

In the limited number of cases Eq. (fTTj) can be solved 
analytically [3 [32|. Numerical methods of solving the 
fractional Fokkcr-Planck equation (fTTj) resume usually 
two directions. First, it is possible to approximate the 
solving probability density p(x, t). This can be achieved 
by numerical approximation to the fractional derivatives 
present in Eq. (JTTj!, see [2(1 E! EH- Nevertheless, due to 
a nonlocal character of fractional derivatives, the conver- 
gence of approximation schemes can be very slow. The 
other possibility relies on the construction of the ensem- 
ble of trajectories from which the estimators of p(x, t) 
can be derived in the form of the frequency histograms. 
Realizations of the stochastic process, whose probability 
density evolves according to Eq. (fTTj) , can be constructed 
either by the subordination technique j3(j, [3l|, Ell Ell or 
by the CTRW framework [37|,[38[. Here we compare these 
two approaches J3l|, [38[ and discuss their rate of con- 
vergence to known, analytical solutions of corresponding 
fractional Fokker-Planck equations. 

In a series of papers [3l|, [39|, 5(| subordination meth- 
ods 0, [30, HJ have been extended to give a stochastic 
representation of trajectories of the process X(t) which 
is otherwise described by the fractional Fokker-Planck 
equation (fTTj. Within the subordination approach the 
process of primary interest X(t) is obtained as a func- 
tion X(t) = X(S v (t)) by randomizing the time "clock" 
of the process X(s) using a different "clock" S v (t) which 
links the real time t with the operational time s. Ac- 
cordingly, the stochastic representation of the solution to 
Eq. (fTTj) is obtained in this scheme by use of the self simi- 
lar process S v (t) representing so called ^-stable (inverse) 
subordinator. The latter is the process with nonncgative 
increments [7| and can be defined as 

S u (t)=m£{s>Q:T(s)>t}. (15) 

In the above equation T(s) denotes a strictly increasing 



instable process (0 < v < 1) whose Laplace transform 
is given by (e~ fcT ( s )) = e~ sk " whereas its inverse S v (t) 
determines random hitting time (first passage time) for 
the problem [y, |7|. In turn, the parent process X(s) is 
composed of increments of symmetric a-stable motion 
64] described in an operational time s by the equation 



dX(s) = dL a , Q (s). (16) 

The combination of Eqs. (fTSj) and (fTTjj) fully determines 
the process X(t) = X(S u (t)) and provides a stochastic 
representation to Eq. Il l I t in terms of the ensemble of 
trajectories [H E! S3 Ha] . 

In a less formal, albeit quite intuitive way, description 
of continuous realization of the CTRW scheme Eqs. Q])- 
(fSj) has been proposed by Fogcdby [29| and Eule [301 ] . 
Their formulation of subordination procedure is based 
on the analysis of the set of coupled Langcvin equations 



x(s) = § = £(,s) 
t(s) = r)(s) 



(17) 



where the random walk x(t) becomes parametrized by 
variable s. In the above equations £(s) and rj(s) are as- 
sumed to be independent, random noises and the pair 
process (x(s),t(s)) preserves the Markov property. The 
requirement of causality (£(s) is a physical time) limits 
choice of ?y(s) to functions returning positive values only. 
The combined process in physical time t is described by 
the trajectories x(t) = x(s(t)) and is subordinated to 
the parent process with corresponding realizations x(s). 
Moreover, the time transformation implies 

s(t) = T <=> t(T) = t, (18) 

or 

Prob {s(t) <T}= Prob {t(T) > t} . (19) 

Let £(£) stands for a white, symmetric Levy noise [44J|45| . 
In this case, the increments Ax(s) are assumed indepen- 
dent over non-overlapping time intervals, i.e. they define 
a stationary process being a generalization of a Brownian 
motion (a generalized Wiener process) . Accordingly, the 
solution of the stochastic differential equation x(s) = £(s) 
can be expressed in the form (NAs' = s) 



x(s) = L a>0 (s) 



r s N-l 
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with the PDF l a fi(x) whose Fourier transform <&(£;, s) 

( e ikx(s)j = ( e ikL a , {s)^ ig givcn by 



$(fc,s) =cxp[-s|fc| Q ] 



(21) 



In turn, let us assume that the stationary random process 
r](s) is a white ^-stable Levy noise which takes positive 
values only. In consequence, the integrated process 



t(s) = / ds'rj(s') 
Jo 



(22) 



is a ^-stable totally skewed Levy motion with an index 
< v < 1 and characteristic function given by 



$(fc, s) = exp 



|fc|" 1 — isignk tan — J 



(23) 



The PDF of the random variable s at time t, p(s, t) is 
then given by the inverse stable distribution 



p(s,t) = -j- / ®{y,s)dy 



{l-L^it/s^)), 
as 



(24) 



with the process s(t) being an asymptotic (continuous) 
analog of the number of steps n(t), cf. Eqs. ©-(JSJ). In 
the above ®(y, s) stands for the PDF of times t with the 
Fourier transform (e~ kt ^) given by Eq. (|23[) . Again, the 
PDF p(x,t) of the subordinated process x(t) coincides 
[29l , |31| with the solution to the fractional Fokkcr-Planck 
equation (fTTj) . 

Let us stress that the numerical methods to approx- 
imate the solution of the Langevin equations (|17p in- 
volve an Ito integral with respect to the (generalized) 
Brownian motion and assume the discretization of the 
time parameter. The computer algorithm generates then 
the approximation of the trajectory (a single realization 
of the process X(t)) in terms of a random walk on a 
one-dimensional grid, i.e. it simulates sample paths of 
a corresponding CTRW. Therefore both, subordination 
and CTRW methods, are essentially different facets of 
the random walk methodology. However, the method 
based on the subordinated Langevin equation provides 
exact representation of solutions of the fractional Fokkcr- 
Planck equation, while CTRW scenarios reconstructs so- 
lutions asymptotically. 

The fractional time derivative in Eq. (fTTj) results in the 
Mittag-Leffler decay of temporal eigensolutions. Conse- 
quently, the alternative approximation [38[ to Eq. (fTT| 
relies on the generation of waiting times distributed ac- 
cording to the complementary distribution function given 
by the Mittag-Leffler function [46[ while the jump lengths 
are distributed according to the a-stable density. 

The probability density p(x,tY see Eq. (fTT)) . is known 
to have a series representation [8[ : 
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where y = x/t v ' a . The above series are divergent for 
a ^ v. However, for a = z/, the summation has been 
shown [8| to produce the closed analytical formula 



p(x,t) 



1 



sm(jrv/2) 



K\y\t \y\» + \y\-» + 2cos(Tns/2) 



(26) 



where (as previously) y = xjt v ' a . In the limit of a = 2 
and v = 1 the standard Markovian (Brownian) diffusion 
is recovered. 



Exact solutions to Eq. (fill) are known in special cases 
f8J, [32|: For v = 1 with any value of a, p(x,t) is an in- 
stable process 



p(x,t) =l a {t 1/a x). 



(27) 



In particular for a = 2 the Gaussian distribution is ob- 
tained 



p(x,t) 
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while for a = 1 the Cauchy distribution is reached 

, , t 1 
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(29) 
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Analytical solutions are also known for v — 1/2, a = 1 



p(x,t) = — 



1 



exp 



2^l 2 ^t 
and for v = 2/3, a = 2 (see Q) 
32/3 
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where Ei(cc) and Ai(x) are the integral exponential func- 
tion and the Airy function respectively. 



III. RESULTS 

Left panels of Figs. [TJ - [6] present sample probabil- 
ity density functions p(x, t) estimated at various times 
t = {1,5,20}. For each parameter set, numerical results 
were constructed separately by the subordination tech- 
nique [U and the CTRW method [H]. Furthermore, 
to ascertain the correctness of both methods, verifica- 
tion of obtained solutions to the FFPE has been exam- 
ined against analytical formula for several known cases 
Eqs. p6] ) . (J28]) . (|2% |l . pp jl and (j3"T]t . Derived results were 
averaged over N = 10 6 realizations of the correspond- 
ing stochastic process. In the subordination method the 
time step of the integration was adjusted to At = 10~ 2 . 
Examples of the histograms generated according to both 
procedures, see Figs. [JJ-[S1 demonstrate that the subor- 
dination method not only reconstructs properly the tails 
of probability densities but also matches correctly their 
central parts. 

The comparison between numerical solutions and their 
analytical analogues (see left panels of Figs. [TJ [6j) is 
based on the analysis of the sum of squared deviations: 



M 1 



N b 



t)f 



(32) 



i=i 



In Eq. (|32|) . Nb represents number of histograms bins, 
Xi locations of bin centers, p(x, t) and p(x, t) denote an- 
alytical and estimated probability density functions, re- 
spectively. Various curves in bottom panels of Figs. [TJ 
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FIG. 1: Probability density functions p(x,t) (left panel) and 
complementary cumulative density functions F c (x) = 1 — F(x) 
(right panel) for the subdiffusion parameter v = 1.0 and the 
stability index a = 2.0. Results were constructed using the 
subordination method (black circles) and the CTRW frame- 
work with the Mittag-Leffler waiting time distribution (empty 
triangles) . Thin solid black line represents analytical solution 
given by Eqs. (J28J) and (|33[) . Various panels correspond to 
various times t — {1,5,20} (from top to bottom). Finally, 
bottom panel presents the sum of squared differences, M 2 , 
between analytical and numerical results (see Eq. Q32[ l) for 
v = 1.0 a = 2.0. Various curves compare simulation results 
constructed by the subordination and CTRW methods with 
various time steps of the integration ( Ai) and number of rep- 
etitions (TV). 



[S] correspond to different method of construction of solu- 
tions or different simulation parameters. 

In general, results based on the subordination tech- 
nique reconstruct much better the shape of theoretical 
probability density. For the Markovian case v = 1, 
and after sufficiently long time, similar accuracy of both 
methods is observed, see bottom panels of Figs. Q] and [5] 
Furthermore, in the Markovian case the level of agree- 
ment between analytical and numerical results depends 



only on the number of repetitions, i.e. both the sub- 
ordination and CTRW methods have similar accuracy, 
see bottom panels of Figs. Q] - [5] On the contrary to 
the Markovian (v = 1) case, for a non-Markov process 
(v < 1), differences between both methods are well vis- 
ible even for long times, see Figs. 02 02 02 and 02 Here, 
the level of agreement depends both on the method ap- 
plied and the number of repetitions. The method based 
on the CTRW with the one-parameter Mittag-Leffler dis- 
tribution of waiting times (with 10 5 or 10 6 repetitions) 
leads to the same level of agreement. In contrast, the ap- 
proach based on the subordination results in significantly 
smaller deviations from theoretical distributions. Also, 
expanding the ensemble of simulated trajectories (by in- 
creasing the number of repetitions) clearly increases the 
level of agreement between analytical and numerical so- 
lutions. In the force free case, which is studied here, the 
choice of the times step of integration seems to be less 
important than the choice of the number of repetitions. 
Finally, the increase in the histogram range (with fixed 
number of bins) leads to smaller values of M 2 because 
spatial resolution of the histogram is decreased (results 
not shown). Nevertheless, such a comparison still demon- 
strates better performance of the subordination method 
than the CTRW framework. 

Right panels of Figs. Q]- El present sample complemen- 
tary cumulative distributions F c (x,t) for various times 
£ = {1,5,20} 

F c (x,t) = l-F(x,t) = 1- / p(x',t)dx' (33) 



with their analytical counterparts obtained after integra- 
tion of corresponding PDFs p(x, t) given by Eqs. (|26|) . 
([55]). (US]), (001) and (EH]). Right panels of Figs. CO -[5] 
demonstrate that both methods perfectly reconstruct the 
asymptotic dependence of the probability densities. 

The properties of solutions to the FFPE (fTTj) are de- 
termined by the subdiffusion parameter v and the sta- 
bility index a. In the simulations based on a subordina- 
tion scheme, the subordinator process S v (t) is evaluated 
by generating first increments of the process T(s), cf. 
Eq. (|T5l) . for which S v (t) forms the inverse. Since T(s) is 
a Levy jump process with nonncgative increments, every 
jump of T(s) can be associated with a long waiting time 
of its inverse S v (t) 0, |47| . This heavy-tailed distribution 
of waiting times is a feature of subdiffusivc dynamics and 
is responsible for a weak ergodicity breaking and a non- 
Markov character of the combined process X(t) [48l - [54| . 

The power law distribution of waiting times (for u < 1 
the mean waiting time is divergent) is also properly re- 
constructed by an explicit use of the Mittag-Leffler dis- 
tribution of the jump-times in the CTRW scheme. In 
fact, the Mittag-Leffler distribution interpolates between 
the short-time stretched-exponential distribution of wait- 
ing times and the long time power-law asymptotics J55| . 
The resulting long-range memory of simulated PDFs is 
well visible in histograms where (for v < 1) a persistent 
cusp at x — is detected (see left panels of Figs. 02 - 



[6]). This behavior typical for subdiffusive systems [56[ is 
also manifested in ambivalent processes like "paradoxical 
diffusion" [43, [5J, [53] ■ In these cases, due to the competi- 
tion between long waiting times (y < 1) and Levy flights 
(a < 2), the second moment of the process X{t) scales 
like (x 2 (t)) oc i 2l/ / Q , i.e. for 2z^ = a it assumes the form 
characteristic for a "normal diffusion", although X(t) is 
non-Gaussian and non-Markov in nature. 

The presence of Levy flights is visible in tails of the 
probability density functions. For a < 2, the comple- 
mentary cumulative distributions F c (x,t) demonstrate a 
power-law decay of the same type like a-stablc Levy den- 
sities governing distributions of jumps (see right panels 
of Figs. El- E]). 

In overall, numerical simulations corroborate that 
asymptotic (space) dependence of the process is deter- 
mined by the jump length distribution while the rate of 
convergence to the long time asymptotics is determined 
by the subdiffusion parameter v, see bottom panels of 
Figs. [T] [6j This is particularly pronounced in the rate 
of decay of differences between theoretical and estimated 
probability densities (M 2 see Eq. (j3"2")l ) when for smaller 
value of the subdiffusion parameter the rate of conver- 
gence is significantly slower. 



IV. NUMERICAL ISSUES 

The continuous time random walk with a one- 
parameter Mittag-Lcfficr distribution of waiting times 
and subordination method provide interesting and effi- 
cient frameworks for construction of solutions to the frac- 
tional Fokker-Planck equation (fTTj) . Despite the fact that 
both methods can be used for solving the same fractional 
Fokker-Planck equation, there are some inherent differ- 
ences among them which we want to discuss in more de- 
tails. 
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FIG. 2: The same as in Fig. [T] for the subdiffusion parameter 
v = 1.0 and the stability index a = 1.0 and Eq. (|29|l . 



A. Precision 

Precision of the constructed results depends on the 
number of repetitions (MLF and subordination meth- 
ods) and the time step of the integration (subordination 
method). Consequently, the method based on the subor- 
dination is more controllable, sec bottom panels of Figs. Q] 

-El 

Figs. [T] - El compare analytical and numerical results. 
These figures indicate that the subordination method 
leads to higher level of agreement between approximate 
and exact solutions p(x, t) of a corresponding fractional 
Fokker-Planck equation. First, subordination method re- 
constructs well not only the asymptotics (tails) of the so- 
lution but also its central part. Next, this method recon- 
structs the analytical solutions for all considered values of 
time t while the CTRW scheme reproduces correctly the 
PDFs only in the asymptotic limits (i.e. after sufficiently 
long times and for sufficiently large space excursions x). 



For v — > 1 the process X(t) becomes a Markov Levy flight 
(see Fig. [2]). Its CTRW approximation scheme is then 
composed by use of the Mittag-Lcfficr function which for 
v = 1 becomes an exponential distribution of waiting 
times. Also in this case, the convergence of both meth- 
ods in reproducing PDFs p(x, t) is met asymptotically. 

Finally, Fig. [7] quantify accuracy of both methods by 
comparing the ratio R of sums of squared differences 
(M 2 ) for the subordination- and the CTRW- formalisms 
at time instant t = 300. Direct analysis clearly indicates 
that the algorithm in which every path is generated as 
a subordination of two trajectories of the processes X(s) 
and S v (t) leads to a higher precision in simulating ade- 
quate solutions to the fractional Fokker-Planck equation. 
The advantage of using the subordination technique is 
especially well visible for small values of the subdiffusion 
parameter v. 
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FIG. 3: The same as in Fig. [T] for the subdiffusion parameter 
v = 0.9 and the stability index a = 0.9 and Eq. (I26I) . 



FIG. 4: The same as in Fig. [T] for the subdiffusion parameter 
v = 2/3 and the stability index a = 2.0 and Eq. (|31|) . 
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TABLE I: Simulation time (in seconds) for various subdiffu- 
sion parameters v and stability indices a. MLF - relates to 
the CTRW framework with the Mittag-Leffler waiting time 
distribution, subordination - to a subordination technique 
discussed in Section [TTJ Results were averaged over TV re- 
alizations with the time step of the integration At. Number 
of repetitions is indicated in the 2nd row while the time step 
of the integration in the 3rd row. 



B. Computational cost 

The simulation time depends both on the number of 
repetition and the time step of the integration. The 
method based on the CTRW with waiting times dis- 
tributed according to the Mittag-Lcmer function is signif- 
icantly faster but less precise. The performance of both 
methods is compared in Tab. |TJ Exact values of the sim- 
ulation time are provided for the informative purposes. 
The ratio between reported simulation times indicates 
which one of the methods is faster. 



C. Extensibility 

The subordination method is easily expendable to 
more general schemes of anomalous diffusion in exter- 
nal fields. In particular, the approach can be efficiently 
used to construct solutions to more general forms of frac- 
tional Fokkcr-Planck equations [40, [58[ with space and 
time dependent forces. On the other hand, the method 
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FIG. 5: The same as in Fig. [T] for the subdiffusion parameter FIG. 6: The same as in Fig. [T] for the subdiffusion parameter 
v = 0.7 and the stability index a = 0.7 and Eq. (|26p . v = 0.5 and the stability index a = 1.0 and Eq. (|30[) . 



based on the CTRW schemes is preferential when model- 
ing anomalous diffusion as stemming from the underlying 
discretized version of a birth-and-death process [59| . 

Both methods can be easily extended to multidimen- 
sional cases. In such realms, conclusions drawn from one 
dimensional case still hold (numerical results not shown). 
However, some additional precaution is necessary when 
generating trajectories of multidimensional Levy flight or 
multidimensional continuous time random walks, as dis- 
cussed elsewhere 1601 • 



V. CONCLUSIONS 

Due to the asymptotic equivalence between various dif- 
fusion schemes and continuous time random walks the 
fractional Fokker-Planck equation plays a special role in 
statistical physics. The continuous time random walk 
scenario with power-law distributed waiting times and 
jumps lengths is asymptotically described by the frac- 
tional Fokker-Planck equation. In majority of situations 




FIG. 7: Ratio R between sums of squared differences M 2 , 
see Eq. ()32[) . for the subordination- and the CTRW- meth- 
ods analyzed for different subdiffusion parameters v. Various 
symbols refer to different exponents a. Data analyzed for 
N = 10 trajectories at time t = 300. 



it is not possible to solve this equation analytically. Very 
efficient and robust numerical approximations to its solu- 
tions are based on the stochastic representation of corre- 
sponding stochastic differential equations driven by sta- 
ble noises (mm. 

In the paper we have examined asymptotic one di- 
mensional diffusion stemming from two alternative ap- 
proaches to CTRW. The method based on the continuous 
time random walks with a one-parameter Mittag-Leffler 
distribution of waiting times provides an efficient way of 
construction of asymptotic solutions. Alternatively, the 
subordination method establishes tools which allow to re- 
produce PDFs of the fractional Fokker-Planck equation 
within the whole time and space domains. Documented 
numerical accuracy and flexibility of the subordination 
method [241 l3l[ 1471 . |63[ results however in higher compu- 
tational costs. 

Altogether, the subordination technique provides a 
useful simulation tool [39| and clearly approximates ef- 
ficiently single time PDFs which solve the fractional 
Fokker-Planck equation. In contrast, the CTRW scheme 
with the Mittag-Leffler waiting time distribution repro- 



duces correctly those PDFs only in the asymptotic limit. 
Our analysis as presented in the paper relates to a 
free diffusion case, when the motion of particles can be 
described by the FFPE (jlip . or alternatively rephrased 
in terms of the CTRW (or Langevin) description (Sec- 
tion |Tl|. This force- free case can be generalized to situa- 
tions with inclusion of the inertia effect. In this context, 
however, special care has to be taken when linking a par- 
ticular subordination scheme with a type of the fractional 
Kramers-Fokker-Planck equation [23|, [3CJ, [47| • 
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